R Notebook

# Import packages
library(tidyverse)
library(lubridate)
library(RQuantLib)
library(plotly)

Step 0: Import data and reference indices; remove superfluous variables to create datasets for agency correction

# Import needed indices
setwd('Datasets/Indices')

The working directory was changed to N:/CommDev/Research/Research/Parks_Regional/Use Estimates/2018/Use Estimate Creation/Datasets/Indices inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

park_ids <- read_csv('Park and Agency Codes Index.csv') # For Step 2

Parsed with column specification:
cols(
  Park_id = col_integer(),
  Park_name = col_character(),
  Agency_name = col_character(),
  Agency_id = col_integer()
)

entrance_classifications <- read_csv('2018 Entrance Usage Classification Index with corrections.csv') # For Step 3

Parsed with column specification:
cols(
  Park_trail = col_integer(),
  Entrance = col_integer(),
  Usage_class = col_character(),
  Start_time = col_integer(),
  Day_type = col_character()
)

bus_vehicle_multipliers <- read_csv('Bus and Vehicle Multipliers Index.csv') # For Step 4

Parsed with column specification:
cols(
  Park = col_integer(),
  ppb_multiplier = col_double(),
  ppv_multiplier = col_double()
)

park_agencies_multipliers <- read_csv('2018 Entrance Usage Classifications with Agency Attribution Index.csv') # For Step 3 & 5

Parsed with column specification:
cols(
  Agency_name = col_character(),
  Agency_id = col_integer(),
  Park_name = col_character(),
  Park_trail_id = col_integer(),
  Entrance_id = col_integer(),
  Start_time = col_time(format = ""),
  Day_type = col_character(),
  Usage_class = col_character()
)

Run the following import and export once, prior to requesting Qualtrics accuracy audits at the end of the season

# Import 2018 Qualtrics datasets - see Qualtrics Info document for login information
# Note that for import by code below to work, all datasets must be csvs, saved with original Qualtrics names in 'Datasets' folder (i.e. download as csv from Qualtrics and do not change the name)
setwd('.')
setwd('Datasets/Qualtrics')

The working directory was changed to N:/CommDev/Research/Research/Parks_Regional/Use Estimates/2018/Use Estimate Creation/Datasets/Qualtrics inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

parks_18 <- list.files(pattern = '*.csv')
make_names <- function(name_vector) {
 
  name_vector_df <- as.data.frame(name_vector)
 
 new_names <- name_vector_df %>%
   rename(Agency_name = 1) %>%
   mutate(Agency_name = str_replace_all(Agency_name, c('\\+-+.+?v' = '', '\\+' = '_')))
 new_names_list <- as.list(new_names)
 
 return(new_names_list$Agency_name)
}
list2env(map(set_names(parks_18, make_names(parks_18)), read_csv), envir = .GlobalEnv)

Duplicated column names deduplicated: 'Q18' => 'Q18_1' [28], 'Q22' => 'Q22_1' [29], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q28' => 'Q28_1' [27], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q21' => 'Q21_1' [28], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.
Duplicated column names deduplicated: 'Q17' => 'Q17_1' [22], 'Q21' => 'Q21_1' [28], 'Q26' => 'Q26_1' [32], 'Q27' => 'Q27_1' [33]Parsed with column specification:
cols(
  .default = col_character()
)
See spec(...) for full column specifications.

<environment: R_GlobalEnv>

all_agencies_list <- list(Anoka_County = Anoka_County, Bloomington = Bloomington, Carver_County = Carver_County, Dakota_County = Dakota_County, Minneapolis_Park_and_Recreation_Board = Minneapolis_Park_and_Recreation_Board, Ramsey_County = Ramsey_County, Saint_Paul = Saint_Paul, Scott_County = Scott_County, Three_Rivers_Park_District = Three_Rivers_Park_District, Washington_County = Washington_County)
# Remove unnecessary variables to clean datasets prior to sending to agencies
# 
# variable_paring <- function(df) {
#   df %>%
#     slice(-2) %>%
#     select(18:33)
# }
# 
# parks_pared <- map(all_agencies_list, variable_paring)
# Write out datasets without extra variables to folder on N drive, from which they can be sent for correction
# parks_pared %>%
#   names(.) %>%
#   map(~ write_csv(parks_pared[[.]], paste0("Datasets/To Send for Accuracy Audit/", ., ".csv")))

Step 1: Tidy agencies’ final data (with corrections), bind together

#setwd('.')
#setwd('Datasets/Audited')
#parks_18 <- list.files(pattern = '*.csv')
all_agencies_list <- list(Anoka_County, Bloomington, Carver_County, Dakota_County, Minneapolis_Park_and_Recreation_Board, Ramsey_County, Saint_Paul, Scott_County, Three_Rivers_Park_District, Washington_County)
park_tidy <- function(df) {
  df %>%
    slice(-1:-2) %>%
  rename(Counter_initials = 18,
         Agency_name = 19,
         Park_trail = 20,
         Entrance = 21,
         Month = 22,
         Day = 23,
         Year = 24,
         Timeslot = 25,
         Scheduled_count = 26,
         Scheduled_date = 27,
         Bikers = 28,
         Peds_skaters = 29,
         Vehicles = 30,
         Buses = 31,
         Horse_riders = 32,
         Boaters = 33) %>%
  select(8, 9, 18:33) %>%
  mutate(Start_time = ifelse(Timeslot != '10:00-12:00 PM', str_replace(Timeslot, '\\-.+? ', ''), '10:00AM')) %>%
  separate(Start_time, into = c('Start', 'AM_PM'), sep = -2) %>%
  unite(Start_time, Start, AM_PM, sep = ' ') %>%
  mutate(Zero = 0,
         Start_time_2 = Start_time) %>%
  unite(Start_alt, Zero, Start_time_2, sep = '') %>%
  mutate(Start_time_fmt = ifelse(Start_time != '10:00 AM' & Start_time != '12:00 PM', Start_alt, Start_time)) %>%
  select(-Start_alt)
  
}
agencies_tidied <- map(all_agencies_list, park_tidy)
agencies_bound <- do.call(bind_rows, agencies_tidied)

Step 2: Create a day type variable (weekend/holiday or weekday); create a variable with the total number of days for that day type in the season

Make sure, in this step, to change the year in the first two lines to the pertinent year - this has to be done manually because while usually use estimates are done the year after summer counts are complete, it could be that they are done in the same year as the counts in years to come

# Calculate days in the summer season for the year
memorial_labor <- as.data.frame(RQuantLib::getHolidayList(calendar = 'UnitedStates', from = as_date(ymd('20180501')), to = as_date(ymd('20180930')), includeWeekends = FALSE))
summer_days <- memorial_labor %>%
  rename(Date = 1) %>%
  slice(-2) %>% # remove Fourth of July
  mutate(row_id = row_number()) %>%
  mutate(Holiday = ifelse(row_id == 1, 'Memorial_Day', 'Labor_Day')) %>%
  select(-row_id) %>%
  spread(Holiday, value = Date) %>%
  mutate(Memorial_Day_Saturday = as_date(ymd(Memorial_Day))- days(2)) %>% # Summer season begins the Saturday before Memorial Day
  mutate(Season_period = as.period(interval(start = as_date(ymd(Memorial_Day_Saturday)), end = as_date(ymd(Labor_Day))), unit = 'day')) %>%
  separate(Season_period, into = c('Season_days', "HMS_in_season"), sep = 'd') %>%
  mutate(Half_of_summer = (as.numeric(Season_days)/2)) %>%
  mutate(Midpoint_of_summer = as_date(ymd(Memorial_Day_Saturday)) + days(Half_of_summer))
# Calculate number of weekend/holiday days and number of weekdays in summer season - this 
day_type_counts <- as.data.frame(seq(as.Date(summer_days$Memorial_Day_Saturday), as.Date(summer_days$Labor_Day), "days")) %>%
  rename(Date = 1) %>%
  mutate(Holiday = RQuantLib::isHoliday(calendar = "UnitedStates", dates = Date)) %>%
  mutate(Day_type = ifelse(Holiday == 'TRUE', 'Weekends_holidays', 'Weekdays')) %>%
  group_by(Day_type) %>%
  count() %>%
  spread(Day_type, value = n)
# Assign which half of the summer the day falls in (this is used, in part, ); assign day-type variable; assign number of days of that day-type in the season (this number is used as a multiplier to scale counts up)
step_II <- agencies_bound %>%
  mutate(Year = ifelse(is.na(Year), '2018', Year)) %>% # Some agencies don't enter the year; fill in for them
  mutate(Month_unite = Month,
         Day_unite = Day,
         Year_unite = Year) %>%
  unite(Date, Month_unite, Day_unite, Year_unite, sep = " ") %>%
  mutate(Date = as_date(mdy(Date))) %>%
  mutate(Half_of_summer = ifelse(Date <= summer_days$Midpoint_of_summer, 1, 2)) %>% 
  unite(Date_year, Month, Day, Year, sep = '-') %>%
  mutate(Date_object = as_date(mdy(Date_year))) %>%
  mutate(Holiday = RQuantLib::isHoliday(calendar = "UnitedStates", dates = Date_object)) %>%
  mutate(Day_type = ifelse(Holiday == 'TRUE', 'Weekend_holiday', 'Weekday')) %>%
  mutate(Total_days_by_type = ifelse(Day_type == 'Weekend_holiday', day_type_counts$Weekends_holidays, day_type_counts$Weekdays)) %>%
  select(-Holiday) %>%
  mutate(Park_trail = ifelse(Park_trail == 'Minnesota River Greenway', 'Minnesota River Greenway Regional Trail', Park_trail)) %>%
  mutate(Agency_name = ifelse(Park_trail == 'Battle Creek Regional Park', 'Ramsey County',
                              ifelse(Park_trail == 'Mississippi Gorge Regional Park', 'Minneapolis Park and Recreation Board',
                                     ifelse(Park_trail == 'Hidden Falls-Crosby Farm Regional Park', 'Saint Paul',
                                            ifelse(Park_trail == 'Scott County Regional Trail', 'Scott County', Agency_name))))) %>% # These agency selection for these trails from the Qualtrics form are incorrect
  mutate(Start_time = as.character(Start_time_fmt)) %>%
  select(-Start_time_fmt)

Step 2a: Add numeric Park Identifier to dataset (used in step 3 to add park entrance usage classifications)

step_IIa <- left_join(step_II, park_ids, by = c('Park_trail' = 'Park_name', 'Agency_name' = 'Agency_name'))
# Ensure all data joined properly
step_IIa_unjoined <- step_IIa %>%
  filter(is.na(Park_id) | is.na(Agency_id))

Step 2a check: Check the number of observations in the step_IIa dataset. If the number of observations increased, some of the counts have been duplicated upon join. Use the following code to discern why and how to ameliorate the faux pas.

# Which agency or agencies have duplicated counts?

# step_II_agency_counts <- step_II %>%
#   group_by(Agency_name) %>%
#   count()
# 
# step_IIa_agency_counts <- step_IIa %>%
#   group_by(Agency_name) %>%
#   count()
# 
# step_IIa_compare_agency_counts <- full_join(step_II_agency_counts, step_IIa_agency_counts, by = c('Agency_name'))
# 
# # In this case, only Anoka.  Which Anoka parks or trails were duplicated?
# 
# step_II_park_counts <- step_II %>%
#   filter(Agency_name == 'Anoka County') %>%
#   group_by(Park_trail) %>%
#   count()
# 
# step_IIa_park_counts <- step_IIa %>%
#   filter(Agency_name == 'Anoka County') %>%
#   group_by(Park_trail) %>%
#   count()
# 
# step_IIa_compare_park_counts <- full_join(step_II_park_counts, step_IIa_park_counts, by = c('Park_trail'))
# 
# step_IIa_compare_counts %>%
#   filter(n.x != n.y)
# # In this case, Rum River Central Regional Park is responsible for the duplication.  Rum River CRP was given two park attributions for Anoka County, by mistake.  Remove the '109' attribution from the index.
# 
# step_IIa_unjoined <- step_IIa %>%
#   filter(is.na(Park_id) | is.na(Agency_id)) %>%
#   select(Park_trail, Agency_name) %>%
#   unique()

Step 3: Add Park Entrance Usage Classifications (High, Medium, Low)

# Note that Trout Brook's RT Entrance 4 was fixed manually; if the 2018 Entrance Usage Classifications Index is re-run in the conversion workbook, this will again need to be changed manually (or scripted) according to the document 'Trout Brook RT Entrnace 4 2018 Designations'
entrance_class_join <- entrance_classifications %>%
  mutate(Entrance = as.character(Entrance)) %>%
  mutate(Start_time = format(strptime(substr(as.POSIXct(sprintf("%04.f", as.numeric(Start_time)), format="%H%M"), 12, 16), '%H:%M'), '%I:%M %p'))
step_III <- left_join(step_IIa, entrance_class_join, by = c('Park_id' = 'Park_trail', 'Entrance', 'Start_time', 'Day_type'))
entrance_class_join %>%
  filter(Park_trail == 61 & Entrance == 1)

Step 4: Add persons per vehicle and persons per bus multipliers, defined by the 2016 Visitor Study and 1998-1999 Visitor Study, respectively

step_IV <- left_join(step_III, bus_vehicle_multipliers, by = c('Park_id' = 'Park'))

Code below checks which parks did not have multiplier

step_IV_unjoined <- step_IV %>%
   filter(is.na(ppb_multiplier) | is.na(ppv_multiplier)) %>%
   select(Park_id) %>%
   unique()
## These parks are new and therefore don't have assigned multipliers; however they have zero vehicles or buses entering the park
# Check full number of observations missing data
step_IV_unjoined_full <- step_IV %>%
  filter(is.na(ppb_multiplier) | is.na(ppv_multiplier))
# Check number of buses and vehicles that came into the new parks
step_IV_unjoined_bv <- step_IV %>%
  filter(Park_id == 133 | Park_id == 140) %>%
  select(Vehicles, Buses) %>%
  filter(Vehicles != 0 | Buses != 0)
## Change missing ppv and ppb multipliers to 1
step_IV_tidy <- step_IV %>%
  mutate(ppb_multiplier = ifelse(is.na(ppb_multiplier), 1, ppb_multiplier),
         ppv_multiplier = ifelse(is.na(ppv_multiplier), 1, ppv_multiplier))

Step 5: Calculate expansion factors & attribute to parks

expansion_factors <- entrance_classifications %>%
  group_by(Park_trail, Usage_class, Day_type) %>%
  count() %>%
  rename(Expansion_factor = n)
# Join multipliers to parks to use to calculate total number of visitors in sample
step_V <- left_join(step_IV_tidy, expansion_factors, by = c('Park_id' = 'Park_trail', 'Usage_class' = 'Usage_class', 'Day_type' = 'Day_type'))
# Check that observations are not missing an expansion factor assignment.  No observation should be missing an assignment (i.e. this code should return a tibble with 0 observations).  If observations are missing an assignment, follow the same steps outlined in the missing assignment check of Step 3 of this script.
step_V_unjoined <- step_V %>%
  filter(is.na(Expansion_factor))
# For purposes of preliminary check of numbers, use only those counts that have an expansion factor assigned (un-assigned factors may be because parks are new or because of incorrect data entry)
step_V_not_na <- step_V %>%
  filter(!is.na(Expansion_factor))
expansion_factors %>%
  filter(Park_trail_id == 114)

Error in filter_impl(.data, quo) : 
  Evaluation error: object 'Park_trail_id' not found.

Step 6: Compute expanded total of visitors (total for entire season, based on sample) at entrance and by date and time of day

# Compute expanded total of visitors (total for entire season, based on sample)
step_VI <- step_V_not_na %>%
  mutate(Expanded_total = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters)) * as.numeric(Total_days_by_type) * as.numeric(Expansion_factor))

Step 7: Add to previous three years of data

x2015_2017 <- x2014_2017 %>%
  mutate(Least_recent_year = min(Year)) %>%
  filter(Year != Least_recent_year) %>% # Filter down to three most recent years of data
  select(-Least_recent_year) %>%
  separate(gooddate, into = c('Month', 'Day'), sep = 1) %>%
  mutate(Year_for_date = Year) %>%
  unite(Date, Month, Day, Year_for_date, sep = '-') %>%
  mutate(Date = as_date(mdy(Date))) %>%
  rename(Agency_id = agency,
         Park_id = Park,
         Bikers = Biker,
         Peds_skaters = Pedskat,
         Boaters = Boat,
         Vehicles = Vehicle,
         Buses = Bus,
         Horse_riders = Horse,
         Usage_class = samclass,
         Total_visitors = totalall,
         Day_type = daytype,
         Start_time = orgstart,
         Expansion_factor = expfact,
         ppv_multiplier = new2008ppv,
         ppb_multiplier = newppb,
         Total_days_by_type = `day#s`,
         Half_of_summer = half,
         Expanded_total = totalsum) %>%
  select(-carpeop, -buspeop) %>%
  mutate(Start_time = str_replace(Start_time, "10000", "1000"),
         Start_time = as.numeric(Start_time)) %>%
  mutate(Usage_class = ifelse(Usage_class == 1, 'High',
                              ifelse(Usage_class == 2, 'Medium', 'Low'))) %>%
  mutate(Start_time = format(strptime(substr(as.POSIXct(sprintf("%04.f", as.numeric(Start_time)), format="%H%M"), 12, 16), '%H:%M'), '%I:%M %p')) %>%
  mutate(Day_type = ifelse(Day_type == 1, 'Weekend_holiday', 'Weekday'))
x2018 <- step_VI %>%
  mutate(Year = 2018) %>%
  mutate(Entrance = as.numeric(Entrance),
         Bikers = as.numeric(Bikers),
         Peds_skaters = as.numeric(Peds_skaters),
         Vehicles = as.numeric(Vehicles),
         Buses = as.numeric(Buses),
         Boaters = as.numeric(Boaters))
  
x2015_2018 <- bind_rows(x2015_2017, x2018)

Error in bind_rows_(x, .id) : 
  Column `Horse_riders` can't be converted from integer to character

Step 8: Compute averages by Park, Trail Entrance/Timeslot Classification, and Day type (Weekend/holiday or Weekday)

x2015_2018_avg <- x2015_2018 %>%
  mutate(Park_id = ifelse(Park_id == 121, 116, Park_id)) %>% # Not sure why, but in past years SW RT was coded as 121.  In all documentation, however, it's coded as 116
  unique() %>%
  group_by(Park_id, Usage_class, Day_type) %>%
  mutate(Average_expanded_total = mean(Expanded_total)) %>%
  ungroup() %>%
  select(Park_id, Usage_class, Day_type, Average_expanded_total) %>%
  unique() %>%
  select(-Usage_class, -Day_type) %>%
  group_by(Park_id) %>%
  mutate(Estimate_18 = sum(Average_expanded_total)) %>%
  select(-Average_expanded_total) %>%
  unique() %>%
  ungroup() %>%
  mutate(Park_id = as.numeric(Park_id)) %>%
  arrange(Park_id)

Step 8a: Compare to previous year’s estimates; note that large differences do not inherently mean incorrect estimates

---

Following is an exploration of outliers

# Plot 2018 counts as striplot; use plotly to mouse over specific observations that appear to be outliers
stripchart_18 <- x2018 %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Park_id, Total_visitors, color = Agency_id)) +
  geom_point(alpha = 0.6, size = 2) +
  coord_flip() +
  scale_color_distiller(palette = "RdPu") +
  labs(title = '2018 Raw Counts',
       y = 'Visitors Counted') +
  theme(axis.text.y = element_text(size = 8))
ggplotly(stripchart_18)

# Plot 2018 counts ; use plotly to mouse over specific observations that appear to be outliers
boxplot_18 <- x2018 %>%
  mutate(Total_visitors = (ppv_multiplier*as.numeric(Vehicles) + as.numeric(Bikers) + as.numeric(Peds_skaters) + ppb_multiplier*as.numeric(Buses) + as.numeric(Horse_riders) + as.numeric(Boaters))) %>%
  ggplot(aes(Park_trail, Total_visitors)) +
  geom_boxplot(alpha = 0.6, size = .5) +
  coord_flip() +
  labs(title = '2018 ',
       y = 'Counted Visitors') +
  theme(axis.text.y = element_text(size = 8))
ggplotly(boxplot_18)